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^v^j ' Abstract. We propose a second order finite volume scheme for nonlinear degenerate parabolic equations. 

For some of these models (porous media equation, drift- diffusion system for semiconductors, ...) it has 
been proved that the transient solution converges to a steady-state when time goes to infinity. The present 
scheme preserves steady-states and provides a satisfying long-time behavior. Moreover, it remains valid 
and second-order accurate in space even in the degenerate case. After describing the numerical scheme, 
we present several numerical results which confirm the high-order accuracy in various regime degenerate 
and non degenerate cases and underline the efficiency to preserve the large-time asymptotic. 
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1. Introduction 

In this paper we propose a second order accurate finite volume scheme for solving the following nonlinear, 
possibly degenerate parabolic equation: for u : IR + X i— > M + solution to 

( d t u = div (f(u)VV(x) + Vr(u)) , x e tt, t > 0, 

(1) 

[ u(t = 0, x) = uo(x), 

with f2 C M. d is an open bounded domain or the whole space E rf , u > is a time-dependent density, / is 
a given function and r € C 1 (JR+) is such that r'(u) > and r'(tt) can vanish for certain values of u. 

A large variety of numerical methods have been proposed for the discretization of nonlinear degenerate 
parabolic equations: piecewise linear finite elements [H HH [301 13H I3S] , cell-centered finite volume schemes 
[2"T1 12"21 l2~i] , vertex-centered finite volume schemes [10] , finite difference methods [33] , mixed finite element 
methods [T] , local discontinuous Galerkin finite element methods [44] , combined finite volume-finite element 
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approach [53]. Schemes based on discrete BGK models have been proposed in [2], as well as characteristics- 
based methods considered in [XT] [32] . Other approaches are either based on a suitable splitting technique 
[20] . or based on the maximum principle and on perturbation and regularization |41j . Also high order 
schemes have been developed in [11, 36, 34 , which is a crucial step getting an accurate approximation of 
the transient solution. 

In this paper our aim is to construct a second-order finite volume scheme preserving steady-states in order 
to obtain a satisfying long-time behavior for numerical solutions. Indeed, it has been observed in [16] that 
numerical schemes based on the preservation of steady states for degenerate parabolic problems offer a 
very accurate behavior of the exact solution as time goes to infinity. To our knowledge, only few papers 
investigate this large-time asymptotic of numerical solutions. L. Gosse and G. Toscani proposed in |28j 
a scheme based on a formulation using the pseudo-inverse of the density's repartition function for porous 
media equation and fast-diffusion equation, and analysed the long-time behavior of approximate solutions. 
C. Chainais-Hillairet and F. Filbet studied in [12] a finite volume discretization for nonlinear drift-diffusion 
system and proved that the numerical solution converges to a steady-state when time goes to infinity and 
F. Filbet analysed the long-time behavior of approximate solution for a king of reaction-diffusion model 
[25] . In [5], M. Burger, J. A. Carrillo and M. T. Wolfram proposed a mixed finite element method for 
nonlinear diffusion equations and proved convergence towards the steady-state in case of a nonlinear 
Fokker-Planck equation with uniformly convex potential. Here we propose a general way for designing a 
scheme preserving steady-states and entropy decay for numerous equations which can be written like ([1} . 

Before describing our numerical scheme, let us emphasize that for some models described by equation 
(fTj). the large-time asymptotic has been studied using entropy /entropy-dissipation arguments, which will 
be the starting point of our approach. On the one hand equation ([I]) with linear convection, namely 
f(u) = u, has been analysed by J. A. Carrillo, A. Jiingel, P. A. Markowich, G. Toscani and A. Unterreiter 
in [7] . On the other hand for equation (fT]) with nonlinear convection and linear diffusion a particular case 
has been studied in [9] [43] by J. A. Carrillo, Ph. Laurengot, J. Rosado, F. Salvarani and G. Toscani. 
We will now remind some of the useful results contained in these papers. 

Case of a linear convection. The paper 7 focuses on the long time asymptotic with exponential decay 
rate for 



Equation |2j is supplemented either by a decay condition when \x\ — > oo if f2 — K d or by a zero out-flux 
condition on <9f2 if is bounded. In the following, we assume that r : K + — > K belongs to C 2 (R + ), is 
increasing and verifies r(0) = 0. We define 



(2) d t U = div (uVV(x) + Vr(u)) , x G H, t > 0, 

with initial condition u(t = 0,x) = Uq(x) >0,u £ L 1 (f7) and 



(3) 




(4) 




and assume that h G L 



loc 



([0,oo)). Then 



(5) 




is well-defined, and H'(s) = h(s) for all s > 0. 

To analyze the large-time behavior to ([2]), stationary solutions u eq of ([2]) in f2 are first studied 




By using the definition @ of h, this can be written as 




and if u eq > in O, then one obtains 



V{x) + h(u eq (x)) = C V.tgO, 
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for some C6l. By considering the entropy functional 

(6) E(u) := / (V(x)u(x) + H{u{x))) dx, 

Jn 

a function u eq ' M e i 1 (r2) is an equilibrium solution of ([2]) if and only if it is a minimizer of E in 

C = ju e L 1 ^), J u(x) dx = M\ 

Under some regularity assumptions on V, existence and uniqueness of an equilibrium solution is proved. 
Therefore, the long time behavior is investigated and the exponential decay of the relative entropy 

(7) £{t):=E{u{t))-E{u^ M ) 
is shown, using the exponential decay of the entropy dissipation 

(8) l{t) := = [ u (t,x)|V(V(s) + h(u(t,x)))\ 2 dx. 

dt Jfi 

Finally using a generalized Csiszar-Kullback inequality, it is proved that the solution u(t, x) of ([2]) with 
r(s) = log(s) or r(s) = s rn , m > 0, converges to the equilibrium u eq ' M (x) as t — s- oo at an exponential 
rate. 

Equation ([2]) includes many well-known equations governing physical phenomena as porous media or 
drift-diffusion models for semiconductors. 

Example 1 (the porous media equation). In the case V(x) — \x\ 2 /2 and r(u) — u m , with m > 1, 
equation f2P is the porous media equation, which describes the flow of a gas through a porous interface. 
J. A. Carrillo and G. Toscani have proved in [10 that the unique stationary solution of the porous media 
equation is given by Barenblatt-Pattle type formula 

l/(m-l) 



m — 1 1 
2m 



x 



(9) u eq (x) = ( d - 

where C\ is a constant such that u eq has the same mass as the initial data Uq. Moreover, the convergence 
of the solution u(t,x) of the porous media equation to the Barenblatt-Pattle solution u eq {x) as t — > oo has 
been proved in [10] , using the entropy method. 

Example 2 (the drift-diffusion model for semiconductors). The drift- diffusion model can also be inter- 
preted in the formalism of 0). It is written as 

' dtN — V • (Vr(iV) - NVV) = 0, 



(10) 



d t P - V • (Vr(P) + PVV) = 0, 



AV = N - P - C, 

where the unknowns are N the electron density, P the hole density and V the electrostatic potential, and 
C is the prescribed doping profile. The two continuity equations on the densities N and P correspond to 
10) with r(s) = s 7 the pressure function. These equations are supplemented with initial conditions Nq(x) 
and Po(x) and physically motivated boundary conditions: Dirichlet boundary conditions N , P and V on 
ohmic contacts T D and homogeneous Neumann boundary conditions on insulating boundary segments T N . 
The stationary drift- diffusion system admits a solution [N eq , P eq ,V eq ) (see [37] ). which is unique if in 
addition: 

(a -is h(N eq \ — V eq i = ° N ^ N cq > h(P eq \4-V eq ) = ap ^ P eq > 

1 ' 1 J \ > a N if N eq = ' 1 ' + \>a P if P eq = ' 

holds, and if the Dirichlet boundary conditions satisfy 111]) and the compatibility condition (if N P > 0) 

(12) h(N) + h(P) = a N + a P . 

In this case the thermal equilibrium (N eq , P eq ,V eq ) is defined by 

AV eq = g(a N + V eq ) - g(a P ~ V eq ) - C on n. 

(is) ; 

N eq = g(a N + V eq ) , P eq = g(a P ~V eq ) on n, 
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where g is the generalized inverse of h, namely 

, . _ ( ft -1 (s) if h(0 + )<s<oo, 

U4j 9[s) -\o if s < h(o+). 

In the linear case r(u) = u, it has been proved by H. Gajewski and K. Gartner in 26 that the solution to 
the transient system U0\) converges to the thermal equilibrium state as t — > oo if the boundary conditions 
are in thermal equilibrium. A. Jiingel extends this result to a degenerate model with nonlinear diffusion in 
[31] . In both cases the key-point of the proof is an energy estimate with the control of the energy dissipation. 

Case of a nonlinear convection. In [S1[5J|33], a nonlinear Fokker-Planck type equation modelling the 
relaxation of fermion and boson gases is studied. This equation corresponds to ([T]) with linear diffusion 
and nonlinear convection: 

(15) d t u = div (xu(l + ku) + V«) , x € R d , t> 0, 

with k = 1 in the boson case and k = — 1 in the fermion case. The long-time asymptotic of this model has 
been studied in ID for both cases [5], in any dimension for fermions [5] and in 3D for bosons [33] ■ The 
stationary solution of (fT5|) is given by the Fermi-Dirac (k — —1) and Bose-Einstein (k = 1) distributions: 

(16) u^(x) = , 

/3e 2 — k 

where /3 > is such that u eq has the same mass as the initial data uq. The entropy functional is given by 

(17) E(u):=[ (^-u + ulog(u)-k(l + ku)log(l + ku)] dx, 
and the entropy dissipation is defined by 



(18) l{t) = I utt + km 



kl 2 , 

— — + loe 



1 + ku 



2 

dx. 



Then decay rates towards equilibrium are given in [9] [8] for fermion case in any dimension and for ID 
boson case by relating the entropy and its dissipation. As in the case of a linear diffusion, the key-point 
of the proof is an entropy estimate with the control of its dissipation. 

Concerning 3D boson case, it is proved in [43] that for sufficiently large initial mass, the solution blows 
up in finite time. 

As explained above, it has been proved by entropy /entropy dissipation techniques that the solution to 
fl} converges to a steady-state as time goes to infinity often with an exponential time decay rate. Our 
aim is to propose a numerical scheme considering these problems and for which we can obtain a discrete 
entropy estimate as in the continuous case. In [3j [6] [5] temporal semi-discretizations have been proposed 
and semi-discrete entropy estimates have been proved. However, when the problem is spatially discretized 
a saturation of the entropy and its dissipation may appear, due to the spatial discretization error. This 
emphasizes the importance of considering spatial discretization techniques which preserve the steady-states 
and the entropy dissipation. This point of view has been already adopted in [TC] but both schemes do 
not provide really satisfying results when the equation degenerates. Indeed both schemes degenerate in 
the upwind flux if the diffusion vanishes and then are only first order accurate in space. Thus we propose 
in this paper a finite volume scheme for nonlinear parabolic equations, possibly degenerate. We focus on 
the spatial discretization, with a twofold objective. On the one hand we require preserving steady-states 
in order to obtain a satisfying long-time behavior of the approximate solution. On the other hand the 
scheme proposed remains valid and second order accurate in space even in the degenerate case. The main 
idea of our new scheme is to discretize together the convective and diffusive parts of the equation ([T]) to 
obtain a flux which preserves equilibrium and to use a slope-limiter method to get second-order accuracy 
even in the degenerate case. 

The plan of the paper is as follows. In Section 2, we construct the finite volume scheme. We first focus 
on the case of a linear diffusion @. Then we extend this construction to the general case (JTJ) . In Section 
3 we give some basic properties of the scheme and a semidiscrcte entropy estimate for the case of a linear 
diffusion ([2]). We end in Section 4 by presenting some numerical results. We first verify experimentally 
the second order accuracy in space of our scheme, even in the degenerate case. Then we focus on the 
long-time behavior. The scheme is applied to the physical models introduced above and the numerical 
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results confirm its efficiency to preserve the large-time asymptotics. Finally we propose a test case with 
both nonlinear convection and diffusion. 

2. Presentation of the numerical scheme 

In this section we present our new finite volume scheme for (|T|). For simplicity purposes, we consider the 
problem in one space dimension. It will be straightforward to generalize this construction for Cartesian 
meshes in multidimensional case. 

In a one-dimensional setting, 57 = (a, b) is an interval of K. We consider a mesh for the domain (a, b), which 
is not necessarily uniform i.e. a family of N x control volumes (Ki) i=1 N such that Ki = x t _ i , x i+ 1 

and 

a = xi < x\ < xs < ... < x,-_i < Xi < < ... < xjv_ < x N ,i = b. 

2 2 2 " 2 J "a:T^ 2 



with Xi — - and 

2 



Let us set 

(19) m( K i) = x i+± ~ for 1 < i < N x . 

Let At be the time step. We set t n = nAt. A time discretization of (0,T) is then given by the integer 
value Nt — E(T/At) and by the increasing sequence of {t n )o< n <N T - 
First of all, the initial condition is discretized on each cell Ki by: 

(20) U? = —^—[ u (x)dx, i = l,...L. 

m(,iViJ J Kf 

The finite volume scheme is obtained by integrating the equation ([T]) over each control volume K t and over 
each time step. Concerning the time discretization, we can choose any explicit method (forward Euler, 
Runge-Kutta,...). Since in this paper we are interested in the spatial discretization, we will only consider 
a forward Euler method afterwards. Let us now focus on the spatial discretization. 

We denote by Ui(t) an approximation of the mean value of u over the cell Ki at time t. By integrating 
the equation |T]) on Ki, we obtain the semi-discrete numerical scheme: 

(21) m{K i )j j U i + F i+ i-F i _ h =Q, 

where is an approximation of the flux — [f{u)d x V + d x r(u)\ at the interface x i+ i which remains to 

be defined. 

Case of a linear convection (f(u) = u). To explain our approach we first define the numerical flux 
for equation The main idea is to discretize together the convective and the diffusive parts. To this 
end, we write [ud x V + d x r(u)] as u [d x (V + h(u))], where h is defined by Qj. Then we will consider 
—d x (V + h(u)) as a velocity and denote by A i+ i an approximation of this velocity at the interface x i+ i: 

(22) A i+i = -dV i+h - dh(U) i+i , 

where dV i+ i and dh(U) i+ i are centered approximations of d x V and d x h(u) respectively, namely 

V(x l+1 )~V(x t ) h(U i+1 ) - h(Uj) 

2 d(Xi,x l+ i) 2 d{Xi,x i+1 ) 

Now we apply the standard upwind method and then define our new numerical flux, called fully upwind 
flux, as 

(23) = F(U h U i+ i) = A+ +i Ui - A- +1 U i+u 

where x + = max(0,x) and x~ — max(0, — x). This method is only first-order accurate. To obtain 
second-order accuracy, we replace in (|23j) Ui and Ui+i by U i+ i _ and U i+ i + respectively, which are 
reconstructions of u at the interface defined by: 

U i+ i = Ui + y^iUi+x-Ui), 



(24) 



U i+i , = U i+1 - U(9 i+ i) (U t+2 - U i+1 ) , 
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with 



Ui+i - Ui 

and <j> is a slope- limiter function (setting = gives the classical upwind flux). From now on we will 
consider the second-order fully upwind scheme defined with the Van Leer limiter: 

_9+\9\ 



(25) 



i + !<?!' 



General case. We now consider the general case where both diffusion and convection are nonlinear in 
Following the same idea as above, we write 



f(u)d x V + d x r(u) = d x [V + h(u)) f(u) 



(26) 

where h(u) is such that h'(u)f(u) — r'(u). Then we define the numerical flux as a local Lax-Friedrichs 
flux 



(27) 

where 
(28) 
and 
(29) 



(f(Ui) + f(U i+1 )) 



A i+i 




2 





(U l+1 - Ui) , 



A i+ i = -dV i+h -dh(U) i+ i, 
a i+ i — max (|/'(u)|) over all u between Ui and L^+i. 



As above, we replace Ui and Ui + i in (f2"7| by reconstructions U i+ i _ and U i+ i + defined by (pM|) to obtain 
a second-order scheme. 



We can now summarize our new numerical flux by: 

A, 



(30) 



?i- 



A i+i 


a i- 











,f(u i+h _) + f(u i+h . 
■■ -dv i+i - dh(U) i+h , 

max(|/'(M)|) over all u between Ui and C/j+i, 



U i+ ^_ = U i + -(t>{6 i ){U i+1 -U i ) 



U i+ i+ = u l+1 - -4 (9 l+1 ) (U l+2 - u l+l ) , 



where either a first-order scheme 
(31) 

or a second order scheme 
(32) 



(j>{6) = 0, 



u. 



u. 



Remark 1 (Generalization to multidimensional case). It is straightforward to define the scheme for 
Cartesian meshes in multidimensional case: the ID formula can be used as it is in any of the Cartesian 
directions. However, the construction of the scheme on unstructured meshes is more complicated. More 
precisely, it is easy to define the first order scheme on such grids, but the difficulty is to obtain high-order 
accuracy. As in the one dimensional case, the idea is to replace the first- order flux F(Ui,Uj), where 



Ui, Uj are the constant values on each side of an edge — K t n Kj 



by F(Uij,Uji), where U l} , U jt 



are second-order approximations of the solution on each side of the edge T^. More precisely, we need to 
obtain piecewise linear functions on each triangle instead of piecewise constant functions. For more details 
concerning these questions, see for example [18, 27 and the references therein. 
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3. Properties of the scheme 

3.1. The semi-discrete scheme. In this part, we study the semi-discrete scheme (j2Tj) - (|30]) - (j32j) and 

consider the equation ([2]) on a bounded domain with homogeneous Neumann boundary conditions. We 
assume that r £ C 2 (R + ) is strictly increasing and h defined by ((U) is in L\ ([0, oo)). Then H given by 
((SJ is well-defined, strictly convex and verifies H'(s) = h(s) for all s > 0. 

We denote by {U? ) i=l N an approximation of the equilibrium solution u eq . This approximation verifies 

(33) dh{U eq ) l+ i +dV l+ i = Vi = Q,...,N x , 
and 

(34) Yl Ax * u * 9 = E Ax * u ? =: M 

j=l »=1 

A semi discrete version of the relative entropy £ defined by ([6]) is given by 

(35) £ A (t) := ^/^(HiUi^-H^-hiUf 9 )^)-!^)). 

i=i 

We also introduce the semi discrete version of the entropy dissipation 

(36) :=E A %^-+i mm [U t+h _(t),U t+h+ (t)y 



i=0 



Proposition 1. Assume that the initial data Ui(0) is nonnegative. Then, the finite volume scheme &21\) - 
if?Tj )) - 01)) for equation 0) satisfies 

(i) the preservation of the nonnegativity ofUi(t), 

(ii) the preservation of the equilibrium, 

(iii) the entropy estimate: for < t\ < ti < oo, 

(37) < £ A (t 2 ) + f'l A {t)dt < £ A (t 1 ). 



Proof. To prove the preservation of nonnegativity, we need to check that 
(38) F (u i+it _,U i+it+ ) - F (u i _ h _,U i _ it+ ) <0 

whenever Ui = 0. 

When Ui — 0, we have U < U+i and U < U-i, and then 6i < 0, which gives <f)(8i) = and finally 

u i+h _ = u i _ h+ = u i = o. 

Then we get 

F (u i+h _,U l+h+ ) F (U^U^) = -AT + JJ i+h+ AtJJ^. 
Moreover, U i _i _ is given by 

which is nonnegative since (f)(6) < 2 for all 9. 

On the other hand, we deal with U i+ i + , and get that either 9i + i < 0, then U i+ i + = J7,+i > 0, or we 
have 9i+i > 0, that is U+2 > U+i and since (f)(8) < 29 for all 9 > 0, we get 

u l+h+ > Ui+! - e i+ i (u+2 - u l+ i) = u l+1 - {Ui+i - Ui) = o. 

We conclude that {55]) always holds when U = 0, which gives (i). 

The part (n) is clear by construction: at the equilibrium, we have dh(U) i+ 1 = 0, which is exactly 

^4 i+ 1 = and then F i+ i = 0. 

By definition ([33]) of £a(*) an d since H'(s) = h(s) for all s > 0, we have 

^(*) = £ a,, (M(c/iW) - nun) 
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Using the numerical scheme (|2"T) . we get 

^( t ) = -j2 mum) - Hu?)) - , 

i=l 

and then a discrete integration by parts yields (using the homogeneous Neumann boundary conditions) 
df N * 

-^(*) = E A ^+i (<w*))<+* - *(^ e3 )i+i) ^+*- 

i=0 

Since by (|33|) we have dh (U^ q ) i+ i = —dV i+ i, we obtain 

= -E A ^^ + i(^ + + i^ + i-(*)-A- + i^ + i, + (*) 

i=0 



< -J^Ax iH \A iH \ min(u i+h _(t),U i+h+ (t)) 



i=0 

Finally we get (Hi) by integrating between t\ and £2- D 

3.2. The fully-discrete scheme. In this part we consider the fully-discrete scheme obtained by using 
the forward Euler method. We denote by U" an approximation of the mean value of u over the cell Ki at 
time t n = nAt. The fully-discrete scheme is given by: 

jjn+l _ jjn 

(39) m( Ki ) ' - i + - J?_ , = 0, 



where the numerical flux F i+ i is denned by (|30[) -(|32 j) . 

Proposition 2. For n>0, assume that U™ > for all i = 1, ...,N X . Then under the CFL condition 
(40) At max IVfci+i) - U^) ~ M^f+i) + M^Dl < l™™Axl 

£/ie fully- discrete first-order scheme i30\) - ![3l\) and L39\) for equation (0) preserves the nonnegativity ofUi, 
which means that U™ + > for all i = 1, N^, and the steady-states solution. 



Proof. Using the definition (|39|) - fj 3 1) - fj31 j) of the fully-discrete first-order scheme, we get for all i — 1, N x 

Ax,; 



Thus we deduce that f/f +1 > as soon as — — is necessarily the 

case from (|4"D|) . using the definition of A n 1 . □ 

Remark 2. 77ms result is not surprising since the stability condition for an explicit discretization of a 
parabolic equation requires the time step to be limited by a power two of the space step. 



4. Numerical simulations 

In this section, we present several numerical results performed by using our new fully-upwind flux. In 
all the numerical experiments performed, the time discretization is given by the forward Euler method. 
We first study the spatial order of convergence of the scheme for linear convection in degenerate cases. 
Then we will apply it to the physical models presented in the introduction: the porous media equation, 
the drift-diffusion system for semiconductors and the nonlinear Fokker-Planck equation for bosons and 
fermions. The results underline the efficiency of the scheme to preserve long-time behavior of the solutions. 
Finally we apply the scheme to a fully nonlinear problem: the Buckley-Leverett equation. 
Below we make comparison between the finite volume schemes (|39p defined with the following numerical 
fluxes: 
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(FU1) 



The first-order fully upwind flux, given by 

A; 



A i+h 








2 





{Ui+x - Ui) , 



with A i+ i, a i+ i denned in (|5D|) . 
• The second-order fully upwind flux, given by 



(FU2) 



1 (f(U, 



A i+i 




2 





U 



u. 



(CU) 



• The classical upwind flux, introduced and studied in [2T] . It is valid for linear convection and 
for both linear and nonlinear diffusion. The diffusion term is discretized classically by using a 
two-points flux and the convection term is discretized with the upwind flux. This flux has then 
been used for the drift-diffusion system for semiconductors [T3J HH [H] . It is defined for equation 
© by 

r (E/i+i) - r (Ui) 



-dV; 



i+4 



Ui 



-dV; 



i+i 



U 



i+l 



Ax, 



i+4 



• The Scharfetter-Gummel flux and its extension for nonlinear diffusion. This scheme is 
widely used in the semiconductors framework in the case of a linear diffusion. It has been proposed 
in 29, 42 for the numerical approximation of the ID drift-diffusion model. This scheme preserves 
equilibrium and is second-order accurate |35j . The definition of the Scharfetter-Gummel flux has 
been extended to the case of a nonlinear diffusion in |16j . For equation © this flux is written 



(SGext) 

where 



dr ; 



A:, 



B 



Ax i+ idV i+ i 



dr 



Ui-B 



Ax, , i dV, 



i+i 



dr. 



U 



i+i 



(41) 



B(x) 



x 



with for a, b £ K+, 



(42) 



dr(a, b) 



e x - 1 
= dr(Ui,U i+1 ) 

h(b) - h(a) 
log(6) - log(a) 

+ b" 



for x ^ 0, B(0) = 1, 



if ab > and a^i), 



elsewhere. 



4.1. Order of convergence. In this part, we test the spatial accuracy of the scheme for linear convection 
(/(s) = s). We consider the equation ([2]) in ID on (— 1, 1) x (0,T) for different values of d x V and r. The 
time step is taken equal to At = 1CP 8 to study the order of convergence with respect to the spatial step 
size. The boundary conditions are periodic. An estimation of the error in L 1 norm at time T is given by 



e2Ax 



\uAx{T) - U2Ax(T)\\ L i(ty, 



where uax represents the approximation computed from a mesh of size Ax. The numerical scheme is said 
to be fe-th order if e 2 Ax < CAx k , for all < Ax <C 1. 

Example 1. We first consider a test case with d x V = 1 and with r{s) = s 2 , thus r'(0) = and r'(s) > 
for all s > 0. The initial data is 

uo(x) — 0.5 + 0.5sin(7ra;), a; €(—1,1) 

and the final time T = 0.1. 

In Table [1] we compare the order of convergence in L 1 norm of the Scharfetter-Gummel extended scheme 
(ISGextj) and of our first and second order fully upwind fluxes (|FU1|) - (|FU2|I . Surprisingly it appears 
that the Scharfetter-Gummel scheme is still second order accurate. This can be explained by the fact 
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that r' vanishes only at one point. Moreover, we verify experimentally that our scheme (|FU2[) is second- 
order accurate and we notice that the L 1 error obtained with it is smaller than that obtained with the 
Scharfcttcr-Gummcl extended scheme. 



N x 


L L error 
SGext 


Order 


L 1 error 
FU1 


Order 


L L error 
FU2 


Order 


100 


1.451.10" 


4 


2 


2.667.10" 


3 


0.87 


8.237.10- 


b 


1.87 


200 


3.619.10" 


5 


2 


1.398.10" 


3 


0.93 


2.208.10- 


5 


1.9 


400 


9.027.10- 


6 


2 


7.156.10" 


4 


0.97 


5.778.10- 


6 


1.93 


800 


2.251.10- 


6 


2 


3.621.10- 


4 


0.98 


1.485.10" 


6 


1.96 


1600 


5.614.10- 


7 


2 


1.822.10- 


4 


0.99 


3.772.10- 


7 


1.98 



Table 1. Example 1 - Experimental spatial order of convergence in L 1 norm. 



Example 2. We still consider equation ^ with d x V = 1, but now with 

r(s) = { (s ~ 1)3 [is ^ h 
^ ' [0 elsewhere, 

then r'(s) = for all s € (0, 1). The initial data is 

Uo(x) = 1 + 0.5 sin(7ra;) x £ (—1, 1), 

and the final time is T = 0.01. 

In Table [5] we compare the order of convergence in L 1 norm of the Scharfetter-Gummel extended scheme 
(jSGextl) and of our first and second order fully upwind fluxes (|FUip - (|FU2p . In this case where r' vanishes 
on a whole interval, it appears that the second-order scheme (|FU2j) is more accurate than the two others 
schemes. The Scharfetter-Gummel extended scheme is only one order accurate while second-order accuracy 
is preserved with our new scheme. 



N x 


L L error 
SGext 


Order 


L 1 error 
FU1 


Order 


L L error 
FU2 


Order 


100 


3.074.10" 


4 


0.96 


2.697.10- 


4 


0.55 


1.053.10" 


4 


1.83 


200 


1.554.10- 


4 


0.98 


1.531.10- 


4 


0.82 


2.830.10- 


5 


1.90 


400 


7.834.10- 


5 


0.99 


8.096.10- 


5 


0.92 


8.040.10- 


6 


1.82 


800 


3.928.10- 


5 


1 


4.163.10" 


5 


0.96 


2.288.10- 


6 


1.81 


1600 


1.966.10- 


5 


1 


2.111.10" 


5 


0.98 


6.576.10" 


7 


1.80 



Table 2. Example 2 - Experimental spatial order of convergence in L norm. 



4.2. The drift-diffusion system for semiconductors. We now consider the drift-diffusion system for 
semiconductors (fT0|) . In the two following examples, the Dirichlet boundary conditions satisfy (fTT |) -([T2 |) . so 
the thermal equilibrium is uniquely defined by (flU)) . We compute an approximation (N!r q , P^ q , V^ q )i = i^..^M x 
of this equilibrium with the finite volume scheme proposed by C. Chainais-Hillairet and F. Filbet in |12) . 

Example 3. Firstly we consider a ID test case on ft — (0, 1). We take r(s) = s 2 . Initial data are 

for x < 0.5 „ „ s f 1 for x < 0.5 



M / \ / fOT X ^ ' 5 P f \ / 

N ^) = { 1 for z>0.5 ' P ^ = { 



for x > 0.5 



and we consider the following Dirichlet boundary conditions 

N(0,t)=0, P(0,t) = l, V(0,t) = -h 
N(l,t) = l, P(l,t) = 0, V(l,t) = l. 

The doping profile is 

-1 for x < 0.5, 



C(x) 



-1 for x > 0.5. 
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The time step is At — 5.10 -5 and the final time T = 10. The domain (0,1) is divided into N x = 64 
uniform cells. 

In Figure [TJ we compare the discrete relative energy £&(t n ) and its dissipation T&{t n ) obtained with the 
Scharfetter-Gummel extended scheme (|SGext[) . the classical upwind scheme (ICUp and our first and second 
order schemes (|FUip - (IFU2[) . The classical upwind flux (|CU[) docs not preserve the thermal equilibrium, 
which explains the phenomenon of saturation observed with it. The Scharfetter-Gummel extended flux 
(|SGext[) preserves the equilibrium at the points where the densities N and P do not vanish, but due 
to the zero boundary conditions on the left for TV and on the right for P, there is also a phenomenon 
of saturation with it. Contrary to these two schemes, our new schemes (|FUl|) - (jFU2|l which preserve 
the equilibrium everywhere, provide a satisfying long-time behavior. Moreover, we computed the relative 
energy and its dissipation with our schemes for different numbers N x of cells and notice that the decay 
rate does not depend on the spatial step size. We obtained satisfying results even for few number of cells. 



Relative energy 



• cu 

* SGext 
o FU1 

□ FU2 



o o o O0 g«Bo n ?«°<' 11 
"n n o n n 



Energy dissipation 



"'i. 



•tittittttiiitttliiitttiiittti 

Q 



• cu 

• SGext 
o FU1 

n FU2 



Figure 1. Example 3 - Evolution of the relative energy E/\(t n ) and its dissipation T&(t n ) 
in log-scale for different schemes (N x = 64). 



Example 4- Let us consider now a 2D test case picked on the paper of C. Chainais-Hillairct, J. G. Liu 
and Y. J. Peng 13 . As in the previous example, the Dirichlet boundary conditions vanish on some part 
of the boundary. The time step is At = 10 -4 , the final time is T — 10 and we compute an approximate 
solution on a 32 x 32 Cartesian grid. 

In Figure [2 we compare the discrete relative energy £A(t n ) and its dissipation lA{t n ) obtained with the 
Scharfetter-Gummel extended scheme (|SGextj) . the classical upwind scheme (|CUj) and the fully upwind 
schemes (|FU1|) - (|FU2[) . We make the same observations as in Example 3: there is a phenomenon of 
saturation with the Scharfetter-Gummel extended and the classical upwind schemes, and not with our 
new scheme. Moreover, the decay rate does not depend on the number of grid cells chosen. 

4.3. The porous media equation. In this part we approximate solutions to the porous media equation 
(43) d t u = V ■ (xu + Wu m ). 

We define an approximation (U^ q )- 1 N of the unique stationary solution u eq (JSJ) by 

where C is such that the discrete mass of {U^ q ) i=l N is equal to that of (Uf) {1 N , namely 

use a fixed point algorithm to compute this constant C . 

i i 
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Relative energy 



Energy dissipation 



B 10 



• 


CU 


* 


SGext 





FU1 


□ 


FU2 



-D?0 o 

""So 
"So 




Figure 2. Example 4 - Evolution of the relative energy E&(t n ) and its dissipation T&(t r 
in log-scale for different schemes. 



Example 5. We consider the following one dimensional test case: m = 5, with initial condition 

J 1 if xe (-3.7, -0.7) U (0.7, 3.7), 
^ ' [0 otherwise. 

Then we compute the approximate solution on (—5.5,5.5), which is divided into N x = 160 uniform cells. 
The time step is fixed to At = 10 -4 and the final time is T = 10. 

In Figure [3] we compare the discrete relative entropy £A(t n ) and its dissipation I&(t n ) obtained with 
the Scharfctter-Gummel extended scheme, the classical upwind scheme and the first and second order 
fully upwind schemes. We obtain almost the same behavior for the Scharfetter-Gummel scheme and the 
fully upwind schemes. We only notice that the dissipation Za(£") obtained with the Scharfetter-Gummel 
scheme saturates before those obtained with the fully upwind schemes. This phenomenon of saturation is 
still greater for the classical upwind scheme. 

Moreover, we compute the discrete L 1 norm of U n — U eq obtained with our second-order scheme. According 
to the paper of J. A. Carrillo and G. Toscani [10] . there exists a constant C > such that, in this case, 

\\u(t) - u e i\\ L1{m < Cexp (-3t/7) , t > 0. 

The experimental decay of U n towards the steady state U eq is exponential, at a rate of about 6, which is 
better than 3/7. 



Relative entropy 





• 


CU 


a 


+ 


SGext 


8 





FU1 


s 


D 


FU2 



O O n * * 



4 6 
t 



Entropy dissipation 



• 


CU 


+ 


SG 





FU1 


□ 


FU2 



Blaesooaiaais 



4 6 
t 



Figure 3. Example 5 - Evolution of the relative entropy £&(t n ) and its dissipation T&(t n ) 
in log-scale for different schemes. 
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Example 6. We still consider the porous media equation, but now in two space dimension on 
(—10, 10) x (—10, 10). We take m = 4 and the initial condition is 



uo(x,y) = < 



exp 

exp 




(- 6 - (a - 2) L ( „ +2)a ) if (x 2) 2 + (y + 2f < 6, 
(- 6 - (a+2) L ( „_ 2)a ) if (x + 2) 2 + (y- If < 6, 



otherwise. 



We compute the approximate solution on a 200 x 200 Cartesian grid, with At = 10~ 4 and T = 10. 

In Figure 0] we compare the discrete relative entropy E&(t n ) and its dissipation I&(t n ) obtained with the 

Scharfcttcr-Gummel scheme, the classical upwind scheme and the fully upwind schemes. 

Figure [5] presents the evolution of the density of gas u computed with our second-order scheme at four 

different times t = 0, t = 0.5, t — 1 and t = 10 and the approximation of the stationary solution u eq 

corresponding to this initial data. 

Moreover, according to the paper of J. Carrillo and G. Toscani [TU], there exists a constant C > such 
that, in this case, 

IKiJ-u^lUiC*) < ^exp(-4t/7t), f >0. 
We compute the discrete L 1 norm of U n — U eq and obtained an exponential decay rate of 2. 



Relative entropy 



• 


CU 


+ 


SGext 





FUI 


□ 


FU2 



*o + 

Bq 4 

H B„ 
no 



Entropy dissipation 



fl « .*** 

am* * 

B 8o 



• cu 



a So„ * 

□ U Q + 

n o • 
□ °o + 

• SGext n n 

"n 5 

o FUI n 

□ 

n FU2 n n 



Figure 4. Example 6 - Evolution of the relative entropy £ a(£™) and its dissipation T&(t n ) 
in log-scale for different schemes. 

4.4. Nonlinear Fokker-Planck equations for fermions. We consider now the nonlinear Fokker- 
Planck equation (|15[) for fermions (k = —1). As in the porous media equation case, we define an ap- 
proximation (U^ q ) i=1 N of the unique stationary solution u eq (|16|) by 

U eq - i = 1 M 

/3e^~ + 1 

where /3 > is such that the discrete mass of (U^ q ) i=1 N is equal to that of ([/?)._ 1 N . We use a 
fixed point algorithm to compute this constant /3. 

Example 7. We consider a 3D test case. The initial condition is chosen as the sum of four Gaussian 
distributions: 

= 27^ l P [ 2— J + ° XP { 2— J + 6XP { 2— J + 6XP v 2— J J ' 

where x x = (2, 2, 2), x 2 = (-2, -2, -2), x 3 = (2, -2, 2) and x 4 = (-2, 2, -2). 

We consider a 40 x 40 x 40 Cartesian grid of fi = (-8, 8) 3 , At = 10~ 4 and T = 10. 

Evolution of the discrete relative entropy £&{t n ), its dissipation I\{t n ) and \\U n — U eq \\^i obtained with 
the scheme (|FU2p is presented in Figure [S] We observe exponential decay rate of these quantities, which 
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(e) Stationary solution 

Figure 5. Example 6 - Evolution of the density of gas u and corresponding stationary 
solution u eq . 

is in agreement with the result proved by J. A. Carrillo, Ph. Laurengot and J. Rosado in [8]. 

In Figure[7]we report the evolution of the level set of the distribution function u(t, x, y, z) = 0.1 at different 

times and the level set of the corresponding equilibrium solution u eq (x, y,z) =0.1. 

4.5. The Buckley-Leverett equation. Finally we consider the Buckley-Leverett equation, with both 
nonlinear convection and diffusion: 

(44) d t u + d x f(u) = ed x (u(u)d x u) . 

The Buckley-Leverett equation is a simple model for displacement of oil by water in oil reservoirs. The 
function u(t, x) represents the fraction of fluid corresponding to oil. The fractional flow function / has a 
s-shapcd form 

u 2 

= 2 i n \2> 

+ (1 — u) z 
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Figure 6. Example 7 - Evolution of the relative entropy £a(£"), the dissipation lA(t n ) 
and the L 1 norm \\U n - U^h. 



and the capillary diffusion coefficient is given by 

v(u) = 4u(I — u). 

The scaling parameter e > in front of the capillary diffusion is usually small. 

Example 8. We consider the following test case [34j|36]: the domain f2 is (0, f), the initial condition 

1-ix if 0<x<±, 



uq(x) 







if | < x < 1, 



and the boundary condition m(0, t) = 1. 

The flux of equation (|44|) can be written under the form (|26|) by taking V — x and 



h[u) = 4 |^log(u) - 3u + 2iT - -u c 

The domain is divided into N x = 100 cells and the time step is At = 10~ 4 . The numerical solution 
computed at different times for different values of e is shown in Figured! The results compare well with 
those in [34l [36] . Moreover, our scheme remains valid for all values of e, even e = 0. In this case the fully 
upwind flux degenerates into the well-known local Lax-Friedrichs flux. 



5. Conclusion 

In this article we have presented how to build a new finite volume scheme for nonlinear degenerate 
parabolic equations. To this end, we rewrite the equation in the form of a convection equation, by taking 
the convective and diffusive parts into account together. Then we apply either the upwind method in the 
linear case or the local Lax-Friedrichs method in the nonlinear case. 

On the one hand, this construction ensures that a particular type of steady-state is preserved. We obtain 
directly a semi-discrete entropy estimate, which is the first step to prove the large-time behavior of the 
numerical solution. On the other hand, we use a slope-limiter method to get second-order accuracy even 
in the degenerate case. 

Numerical examples demonstrate high-order accuracy of the scheme. Moreover we have applied it to 
some of the physical models for which the long-time behavior has been studied: the porous media equa- 
tion, the drift-diffusion system for semiconductors, the nonlinear Fokker-Planck equation for bosons and 
fermions. We obtain the convergence of the approximate solution to an approximation of the equilibrium 
state at an exponential rate. A future work would be to prove this exponential rate by using a discrete 
entropy /entropy dissipation estimate as in the continuous case compared with previous approaches. 
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(e) t = 10 (f) Stationary solution 

Figure 7. Example 7 - Evolution of the level set u(t,x,y, z) = 0.1 and level set of the 
corresponding stationary solution u eq (x,y, z) = 0.1. 
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